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ABSTRACT 


The effects of point defect Pent ae ion in copper and tungsten 
crystal lattice have been studied by computer simulation techniques. 
Vacancies, interstitials, and replacement impurities have been 
created in the first five layers of the free (100) surface of these 
crystals. The subsequent binges enercies of these defects in 
tungsten were compared with experimental temperature dependent de- 
sorbtion peaks, corresponding to binding energies of neon defects 
in a tungsten crystal. tea eet eee and replacement impurity 
positions in the first three to five layers were found that seem 
to correspord to the experimental data. Significant results were 
also obtained which were associated with general surface effects, 


especially crowdion migration. 
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I. INTRODUCTION 


Extensive research has taken place in the last decade in the 
area of computer simulation of radiation damage in crystal lat- 
tices. Two major areas of simulation have been defined. "Dynamic 
Simulation" suggests the firing of an atom or ion against a crystal 
and the observation of the resulting many-body collisions. Ex- 
amples include sputtering simulation Pe 23a in which atoms are 
ejected from the surface of an 1on-bombarded crystal; and chan- 
neling simulation 4], in which ions are fired down open channels 
in Menecllose-packed structures such as body-centered cubic and 
moamend lattices. "Static simulation", on the other hand, is con- 
cerned with equilibrium positions and energies for point defects 
in crystals bee or ale This latter area was the concern of this 
Wee car ch. Specifically, this simulation attempted to correlate 
equilibrium potential energies of point defects with experimentally 


determined binding energies of point defects in Tungsten Roo 14 


A. HISTORICAL BACKGROUND 
1. The Problem 

Historically, all Bee tal dynamics computer simulation has 
been based on the assumption that the complicated many-body problem 
can be reduced to many two-body problems. This assumption has re- 
peatedly been shown to be a valid one when employed incrementally. 
Incremental calculations are necessary Since a complete solution 
in closed form is impossible. Small time increments At approximated 


the true time differential dt of the impossible closed form 





Solution. Specifically, the desired type, size, and orientation 
of crystal lattice is stored in the comouter, appropriate inter- 
atomic potentials are chosen, and all mutual forces between atoms 
are calculated, based on the analytic potential functions. Point 
defects are introduced, and each atom is then allowed to move in- 
crementally, based on these forces and Newtonian Mechanics (5). 
Through proper choice of time increment duration, damping of forces 
and velocities in each time increment, and sufficient repetition 
of the procedure, realistic results are obtained. 
peice PtoOneers and Their Contributions 

Pioneering work in this field began in the early 1960's 
at the Brookhaven National Laboratory. Gibson, Goland, Milgram 
and Vineyard (GGMV) [5] published the results of extensive work 
PMONStAtHe Gia Gynamic Simulatison. In etetic simulaticn the 
determined equilibrium positions for interstitials and associated 
potential energies of formation. In dynamic simulation they in- 
vestigated momentum propagation directions of energetic knock-on 
atoms (focusing), collision chains, and related topics. They used 
a central-difference method to obtain velocities and positions frem 
calculated forces. All work was done with copper, and results were 
correlated with experimental data. Johnson and Brown (JB) [6] did 
extensive work in static simulation, again with copper. They es- 
tablished that only one stable position exists for a single face- 
centered cubic (FCC) self-interstitial: the (100) split inter- 
stitial. Johnson [10] later published further work in this area, 
with formation and activation energies for various point defects. 


Enginsoy, Vineyard, and Englert (EVE) (7) and Johnson [11! 





repeated most of the earlier calculations in GGMV and Johnson 
for the body-centered cubic (BCC) case, based on @ iron. They, 
too, established the existence of only one stable interstitial 
Pesition, a (110) Spliteemeerstitial. 


Girifalco and Weizer (GF) [12] calculated Morse Potential 


one = pl exp (-20(r 5 -r 5) }- Z expt -a(r s -x 3] 


parameters for various metals, based on experimental values for 
the energy cf vaporization, the lattice constant, and compressi- 
bility. Resulting elastic constants and equations of state agreed 
satisfactorily with experiment. Girifalco and Weizer (13) later 
published results of using these Morse parameters in simulating 
vacancy relaxation dynamics. Anderman C14) used GW's technique 


ee eel 


() 


ulating Morse parameters, but instead of summing over an 
entire crystal, (GW calculated out to the 150th nearest neighbor ) 
Anderman found parameters as a result of summing out to second, 
third, and fourth nearest neighbors, for use in short-range approxi- 
mations. 

Harrison [1,2,3] has investigated sputtering phenomenon and 
other surface effects with a modified Brookhaven model, the most 
Significant change being the use of an average force method {15] 
instead of the central difference method in integrating the 
equations of motion. (See Appendix C.) He has also calculated 
repulsive potentials of the Born-Mayer type SNe = Se 
for many ‘combinations of atoms and ions based on secondary elec- 


tron emission, and Hartree-Fock atomic electron distributions tio). 





Beeihne metential Function Problem 

The most difficult problem encountered by computer sim- 
ulation has been the proper choice of the potential function. No 
Simple analytic expression, based on either theory or experimental 
data, has ever been found that completely describes crystal dy- 
namics fie) although many analytic expressions are partially cor- 
rect. The problem has been three-fold: First, present analytic 
expressions have narrow regions of validity, i.e., some correctly 
describe atomic behavior at equilibrium distances, but fail at 
Shorter or greater distances. Second, some analytic expressions 
are limited because they only apply to interactions between iden- 
tical atoms. ere the assumed functions have spherical symmetry, 
and are technically limited to interactions between closed shell 
atoms or ions [18]. Although our assumption of a spherically sym- 
metric potential in crystals is only approximately correct, it 1s 
nevertheless a very good approximation for FCC structures, anda 
reasonably good approximation for BCC structures. It 1s grossly 
in error when applied to diamond structures. 

The atomic potential, with the familiar potential well, 
sharply repulsive wall and Eee attractive tail, varies greatly 
between different pairs of atoms: Since Eee coae give only ap- 
proximate parameters for this complete potential function, experi- 
mental data have been used extensively in the formulation of po- 
tentials. Other avenues have been opened by computer simulation. 
Since potential well depths are typically on the order of a few 


eV, the characteristics of the well can be ignored in high energy 


dynamic simulation. Low energy dynamic simulation, and even static 


10 





simulation, have also been based upon this approximation with 
useful results. Historically this is how crystal simulation 
began. GGMV [5] employed a purely repulsive potential of the 
Born-Mayer (BM) type and applied external forces on all crystal 
boundaries to hold the crystal together. JB [6] used basically 
the same technique. For improved equilibrium studies a potential 
with a well was necessary so GW (213) used a Morse potential in 
their simulation. The Morse function, however, fails at strongly 
repulsive distances. To satisfy the need for a more versatile 
potential, capable of handling both high energy and near-equili- 
brium dynamics, composite potentials were developed, which re- 


semble BM or Bohr 


t 


ee 
(Vi, = yo exp(AtBr, ,)) 


1) 
functions at short separations, and Morse functions at equili- 
brium and greater separations. Specifically, EVE 17] combined 
a screened Coulomb or Bohr potential, a BM potential, and a Morse 
potential, in the higher repulsive, lower repulsive, and attractive 
regions, respectively, of the atomic potential. 
Johnson [10,11] in his later papers, used three cubic equations to 
approximate the meet entia lL. Anderman [14] and Harrison [1,2, 
3,4] have used the BM repulsive term together with a Morse well 
and attractive tail, smoothly fit together by a cubic equation in 
the region near their intersection. 
4. The Point Defect Problem 
All early simulation was done with homogeneous systems: 
All atoms were exactly the same, limiting high energy dynanic 


simulation to bombardment by atoms identical to the lattice atoms, 


11 





and static simulation to consideration of only the vacancy and 
self-interstitial cases. This limitation was forced by the po- 
tential functions, because parameters for the Morse function were 
based on experimental data for homogeneous media (ae). The 
methods used could not yield parameters for different-atom pairs. 
BM parameters, however, are obtainable for different-atom pairs, 
by methods such as the Hartree-Fock method (16) mentioned Dpie- 
viously. 

In spite of these limitations, dynamic computer simulation 
of bombardment by foreign atoms, or static simulation of foreign 
interstitials, can be done by two alternate methods. First, 
Harrison [4] has neglected the attractive interactions and has 
done foreign particle dynamics using repulsion only. Alternately, 
Johnsen [il] has dorived a cubic equation for a complete potential 
with a potential well, based on limited experimental data on car- 
bon defects in iron. However, experimental substantiation for a 
foreign-particle potential well is much more difficult than for 


an identical atom potential well. 


B. THE EXPERIMENT 
The experimental data which this simulation proposed to explain, 

were published by Kornelsen and Sinka (KS) L8]. They have bonm- 

| . + + + + 
barded a clean (100) tungsten surface with Ne , Ar , Kr , and Xe , 
in the energy range of 40 eV to 5 keV. The subsequent "damaged" 
crystal was heated at a constant rate, and gas desorbtion rates 
were measured. Instead of a constant desorbtion of ions, various 


distinct peaks were found, categorized into two basic types: a 


iz 





Single large peak at 1800°%, the same for all four ions; and four 
or five smaller peaks in the 400 °K to 1650°K range, which were 
not in the same position for all four lions. These latter peaks 
were postulated to correspond to binding energies of various point 
defects in the first few layers of the tungsten crystal. (See 
Figure 4.) 

This simulation used Harrison's assumption of a repulsive 
potential only, for interactions between a foreign point defect 
and other atoms in the lattice. When investigating neon defects 
in tungsten, all tungsten lattice interactions were based on com- 
posite Morse and repulsive Born-Meyer potentials, and all neon- 
tungsten (Ne-W) interactions were based on a purely Penenteaye BM 


potential. 


nS 





ieee OBIECL IVE 


The long-range objective of this simulation was to correlate 
simulated and experimental binding energies of neon point defects 
in tungsten. Since the assumption that all Ne-W interactions are 
purely repulsive was not realistic, the degree to which subsequent 
Simulation results are valid must be based on a known standard. If 
the simulated binding energies are not correct, a valid correction 
factor can be applied, if derivable from the known standard. One 
standard which proved to Pele this information was the tungsten- 
tungsten (W-W) interaction. A tungsten point defect could be 
treated as an atom of the lattice, and given an interatomic po- 
tential identical to all other lattice atoms, the composite Morse 
and BM potential. This was Method 1. A tungsten defect could 
also be treated as foreign, and allowed to mLcraee with other 
lattice atoms with a repulsive potential only (Method 2). If a 
specific tungsten point defect 1s treated by both of these methods, 
an empirical relationship between the repulsive potential assumpt-. 
lon and the "true" potential for W-W interactions is obtained. 

The objectives of this research were fourfold: 

1. Demonstrate that the two methods of treating tungsten point 
defect in a tungstenlattice yield basically the same physical 
results, and agree with published results yey concerning 
split interstitial positions. 

2. Develop a general empirical relationship between the binding 


energies derived by the two methods. 


14 


3. Obtain values for binding energies of neon defects in all 
possible positions in a tungsten surface. Transform these 


values to more realistic ones using the empirical relationship 


gemived in 2. 


4. Compare these results with KS's experimental data. 


1 
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Pit. InE MODEL 


Erm tii: CRYSTAL 

The model used in this research is the Gay-Harrison [19] model, 
with modifications by Levy m2 One Johnson baie Effron aie and 
Moore {23]. Abbreviations in brackets refer to computer program 
names for the variable in question. 

Both copper and tungsten crystals were simulated. Copper was 
Simulated only to provide an interface between this research and 
published simulation results. Copper forms a face-centered cubic 
crystal with an experimentally determined lattice constant, (Ee) 
or cube edge distance of 3.615A. The lattice unit (LU), defined 
as 4LC. is 1.80758: and the nearest neighbor distance, as in all 
FCC structures, is \Y2LU. Tungsten forms a body-centered cubic 
mecal, Wage erase G 2o 1 3.164, a LU of 1.58A, and a nearest neighbor 
G@astance, peculiar to all BCC structures, of V 3LU. All distances 
in the program are measured in LU. The program could construct 
feeoo) (110), and (111) orientations of face-centered and body- 
centered cubic structures. wine copper crystal size was 8 x 8x 8, 
and contained 256 atoms for the (100) orientation. 

The major portion of the simulation was done on the (100) 
Orientation of tungsten, corresponding to KS's experimental work. 
This tungsten crystal size for Neon point defects was 10 x 10 x 10, 
and contained 250 atoms. Some W-W simulation was done ona 
14 x 14 x 14 crystal; the reasons are explained in RESULTS. The 
bottom two layers of the lattice were not allowed to move, al- 


though they had potential energy, and exerted force on all atoms 


16 


in the crystal. The other eight layers were completely free to 
move, and were included in the dynamic calculations of each time- 
step. 

The surface layer (Y = 0) and the second layer (Y = 1) were 
moved forward, simulating actual surface relaxation in the crystal. 
This relaxation was calculated by Moore (24) using simulation 
techniques, and tested against previous results by Burton and 
Jura [25]. Definitions and use of mobile layers, relaxation, etc., 
were analagous in the copper model, as were all other aspects of 


the model to be described in this chapter. 


B. THE POTENTIALS 
1. The W-W Composite Potential 
The attractive potential used was the Morse potential, with 
tungsten parameters calculated by GW f12]. The interaction eneray 


oS ee 


ore Stra pair eor particles 1 and j is: 


955 = DLexp{ -20(r 5 .-r ,)} - 2 expl-a(rs-ro)} J (1) 


where D[DCON] is the dissociation energy of the pair, r URE}. is 
the equilibrium separation, x, jlbrst] is the actual separation, 
and @ {ALPHA | 1s a constant. 

The repulsive potential is of the BM type, with Harrison's 


Hartree-Fock parameters. The interaction energy ee 1s 


Viz = exp (A+Br | 5) (2) 


where B [EXB] is always negative, A Lexa] is always positive, and 
Tay (prist] is the actual separation. The constants [EXA] and 


CExB] are peculiar to the W-W interaction. 
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The ranges of the W-W composite potentials were as follows: 
the BM repulsive potential operated from O to 1.58: and the Morse 
potential from 2A [ROEB] to 5.384 [ROEC]. In LU, the dimension 
in which all calculations were done, these constants were .9494, 
1.2658, and 3.4000 Lu. [ROEC] was chosen to include interactions 
out to the fourth nearest neighbor (NN4) at \f1l LU = 3.317 LU but 
not NN5 interactions, at V12 LU = 3.464 LU. Note, however, that 
Slight displacements of NN5's might allow their inclusion in po- 
tential and force calculations. The gap between 1. 5A and 2A was 
fared with a cubic function, which matched to the other po- 
tentials and slopes at [RozA] and LRogeB]. 

2. Purely Repulsive Potentials 

For foreign point defect interactions, i.e., Ne-W, or 
W-W Method 2. a repulSive potential only was used. The potential 
was again a BM, with the constants labeled ([PEXA] and [PEXB]. For 
the W-W, Method 2 interaction, LPEXA] = [LEXA] and [PExB] = [ExB). 
Ranges for foreign point defect interactions, however, were dif- 
ferent, and the potential itself was modified at the cutoff point. 
Whereas the BM part of the composite potential extended out to 
about .95 LU, the modified BM potential used for foreign defects 
was allowed to extend to Vi bu, corresponding to the NN1 distance 
[RoE]. Cutting the potential off at WROnietert a.step Of about 
.O5 eV for Ne-W (.2 eV for W-W) at the NN1 equilibrium position. 
Since neither discontinuities nor repulsive potentials were de- 
sired at this equilibrium position, we "eroded" [15] the potential 
by subtracting v(CROE)), or about .05 eV from ae Lox 
oe 2 reel. Calculated forces, based on these eroded potentials 


must wbe modified also, but for a different reason. It is possible 


18 





to conceive of a case where an atom is further away from the 
defect than [ROE] at the beginning of a timestep, but closer than 
[ROE] at the end. The force has essentially "turned on" in the 
middle of the timestep. The modification gives such an atom a 
force which is less than the final force by approximately a fac- 
tor proportional to the ratio of distance traveled outside [ROE] 
to the total distance traveled during the timestep a. (See 


Appendix B.) 


See foe TIMESTEP 

Motion caused by these forces must be found by an approximate 
numerical method of time integration. As in previous work with 
this model, the average force method [15] was used. In this 
method, all mutual forces were calculated in subroutine STEP. 


meeeu- om these <Oorces, new tempcrazy velocitizes and nositions 
5 t By fc 


were found. Forces were again calculated, based on the temporary 
positions. The final positions were then calculated from the 
average of these two force determinations. All velocities were 


then either zeroed or halved as a damping method. This consti- 
tuted one timestep. 

For this average force method to work properly, the atl pr} 
which approximates dt in the integration must be kept small. Too 
Small a value for fone however, would result in excessive com- 
puter time. The choice of [pr] was also complicated by the fact 
that Cpr] must be kept much smaller earlier in the program, when 
Velocities, forces, and pretgics are large, but can be allowed to 
grow larger as the simulation approaches equilibrium. For this 


reason, at the end of each timestep, a new [pr] was calculated 
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for use in the next timestep. The parameter chosen to control 
(pr] was (przx], the distance, measured in LU, which the most 
energetic atom was allowed to move before starting a new timestep. 
[prij has varied between .001 and .02, depending on such con- 
ditions as otiginal position of the point defect, relative masses 
of atoms, etc. In general, ( pri) must be kept very small when 
high velocities are expected, and can be increased when all motion 
1s expected to be "sluggish". In actual practice lpmivand prt 
are related to both the velocity of each particle and the force 
on each particle. To insure that no particle traveled more than 
Corr] we ensured that Cpr] was small enough so that neither the 
velocity of the most energetic atom nor the force on the most 
stressed atom would result in motion greater than pra): (See 


Appendix C.) 


D. FOREIGN INTERSTITIALS 

1. Unequal Mass Implications 

Many changes in the program were necessary when a foreign 

defect was included in the lattice. These changes were especially 
necessary when the defect was much lighter than the lattice atoms; 
1.e€., neon in a tungsten lattice. First, in the average force 
calculations, a separate section had to be added for the calcu- 
lations for the primary, or "bullet'', based on the bullet mass 
LpMAS]. Second, the potential energy between two unequal mass 
atoms was split in proportion to their reduced masses (See 
Appendix D). Third, the section that determined the new timestep 
duration was originally based cn the lattice atom mass. Since the 


light interstitial is usually the most energetic or most stressed 
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atom, very erratic behavior was observed until the timestep du- 
ration calculations were revised to handle two different masses 
(see Appendix C). Finally, a significant mass difference between 
defect and lattice atom required a reduction in [pri]. For the 
neon-tungsten simulation, a Lori] of .5% was used. 

2. Tonization State and Repulsive Potentials 

The major portion of the Ne-W work was done with the 

assumption that tungsten was ina +6 state, and neon was neutral 
in the lattice. Experimentally KS fired neon in a +l state into 
tungsten, but once emplanted, the ionization of neon was unknown. 


+] +6, 4 
and W with Neo and N were sub- 


All combinations of Ww, W 
. : 4 + 
jected to Hartree-Fock analysis. (See Figure 5.) Only Ne-=-W 
interacted in an approximately exponential manner and could there- 
fore possess realistic BM parameters. Attemnts to linearize 

© .,+1 O06 
Ne -W and Ne -W were made, and subsequent BM parameters were 


determined. The results of such changes did not Significantly 


influence the results of this investigation. 


E. RUNNING TIME 

The following factors effected the problem Ganda time: range 
of potential, size of crystal, depth of mobile layers, and degree 
of damping. First, the range of the potential was picked to 
include at least NN4 interactions. The range used for the tungsten 
Simulation was 3.4 LU, which includes interactions out to NN4&. The 
error made by neglecting NN5 interactions was only 3% in the 
binding energy of an interstitial, but the omission of NN5 inter- 
eens) cit running time almost 10%. Second, the size of crystal 


and depth of mobile layers were picked as small as possible, fcr 
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reduced running time, but were at least large enough to completely 
contain the potential range. Third, the half velocity method of 
damping was used whenever possible. In general, when velocities 
were zeroed at the end of each timestep, a timestep took about ten 
seconds, and equilibrium was reached in about 300 timesteps. When 
velocities were halved, each timestep again took about ten seconds, 
but equilibrium was reached in about 150 timesteps. 

It was originally expected that increasing [pri] would decrease 
running time. Most W-W simulation was done with Vipsesei) 29.) Vets 
the most energetic or stressed atom could travel .02 LU before the 
damping of velocities and the starting of a new timestep. When 
ort) was increased, the atoms moved more erratically toward equi- 
librium, and vibrated there, but did not achieve equilibrium 


Seeeeetcantly scorer. (See Pisure ¢.} 


F. SUMMARY 
In summary, the steps of the program are outlined: 

1. Variables are initialized, constants established, and input 
data read in. 
Zeocaling factors and time Saving multipliers are calculated. 
3. Morse and BM potential functions are calculated based on 
input data. Subsequent forces, based on derivatives of these 
functions are calculated. Potential erosion and force modifi- 
cations are performed. 
4. Potential cutoff's [RoEA], [ROEB], LROEC] are established 
and the smooth fitting cubic equation is placed in the gap. 
Dis Phe anes ccd crystal type, size, and orientation is built, 


and the point defect positioned (see Appendix A). 


ee 





6. Mutual potential energies of all atoms in the crystal are 
calculated. Local potential energy is calculated. (See 
Appendix D.) 

7. All initial positions and potential energies are printed, 
along with total potential and total kinetic energy, local 
potential energy, and the change in local potential energy. 

8. The first timestep is started, with an arbitrary running 
time of 10714 sec. Velocities and positions are calculated 
by the average force method, and the maximum velocity Cemax] 
and maximu;n force [FMAX] are found. 

9. A new [pt], based on CEMAX], CFMAx], and [pri] is calculated 
for use in Ate next timestep. (See Appendix C.) 
10. All velocities are zeroed or halved as an energy damping 
method, and the process (8. to 10.) is repeated. 

11. At selected timesteps, all changes in position (Cpdx], Coy], 
| (pz), velocities (vec Lvy], [vxl), and kinetic, potential, 
and total energies [PKE], (pPE], [pte] for each atom in the 
crystal are printed. 
12. The program is ended after a pre-selected timestep, with a 
final printout of Bocuiton and potential energy of each aton, 


as in Step 7. 
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Ve ORE SULLS 


KS's experimental data indicated that four or five interstitial 
positions in the first few layers of a tungsten lattice could be 
found that would result in different binding energies. It was soon 
found that many parameters in the program could effect the results, 
and so a systematic attempt to isolate the effects of each indivi- 


dual parameter was undertaken. 


Aa THE CRYSTAL 

As explained in Appendix A, the Y = O plane was the crystal 
surface; the Y = 1 plane was the first layer beneath the surface, 
‘etc. Each atom in each layer was then designated by appropriate 
x and z coordinates. This construction was independent of type 
mamleattice- 2.e€,, the surface layer in either BCC (100) or FCC 
(lOO) "was Y = 0, etc. An interstitial that escaped the lattice 
normal to the surface travelled in a (010) dizectien, and am 26m 
that escaped normal to a side travelled in either a (100) or (001) 
direction. 

The dimensions of the lattice had a great bearing on the re- 
sults. In general, Ehemlatge. the crystal, the more realistic the 
results, but increased computer time prevented the use of a size 
bigger than absolutely necessary. All point defects were placed 
as close to the center of each plane as possible, and the x and z 
dimensions of the lattice [1x] and [1z] were chosen to completely 
enclose a circle of radius [ROEC] from the point defect. In this 


way any point defect in the center of the lattice would not feel 
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the effect of the sides of the lattice, especially unequal numbers 
of atoms in all directions. [In the y direction, the lattice was 
again built deep enough to completely contain the radius of the 
potential of a point defect placed at the center of the lattice. 
To simulate the effect of an infinitely deep lattice, the bottom 
two layers of the crystal were held immobile, but still allowed 
to interact with all mobile atoms above them. The tungsten crystal 
size used most often was a 10 x 10 x 10 cube with the bottom 
items x 10 volume held rigid. 

Often erratic behavior in the simulation could be eliminated 
by Simply increasing the crystal size. This was especially true 


for the problem of crowdion migration. 


B. CROWDION MIGRATION 

Growdion Mfagration 1s a chain reaction of single lattice site 
jumps initiated by interstitial implantation. If the chain re- 
action ends by pushing the surplus atom into an already existing 
vacancy, the interstitial-vacancy pair is called a Frenkel pair. 
A Frenkel pair can also be created dynamically by moving an atom 
fuem 1ts lattice site to a nearby interstitial Roshi on, from 


where it can cause migration back to the vacancy. If the mi- 


gration cannot find a vacancy, and travels all the way to the 


surface, the surplus atom forms a "stub'"'. Normally, migration 
is always in a closed packed direction; i.e., in the Cini): dae 
mection in BCC. It was discovered, however, that this rule wes 


modified near a surface, since an imbalance of forces in the di- 
rection normal to the surface automatically pushed a crowdion in 


that normal direction into a stub position. [n the tungsten 
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lattice, for instance, a tungsten interstitial that did not 
initiate crowdion migration would reach equilibrium in a (110) 

Split interstitial position, as previously found by EVE eek 
tungsten interstitial that did initiate crowdion migration would 
sometimes migrate in a (111) direction because of closed-packedness , 
or sometimes in a (100) direction if implanted near a (100) sur- 
face. Crowdion migration was never found in a (110) direction, 
Since this is the least closed-packed of these three directions. 
(Hence the tendency toward split-interstitials in this 

direction). 

Crowdion migration was a very common process near a lattice 
surface. It was found, however, that varying the ae of atoms, 
the range of the potential, and the rate of energy damping could 
Pere Col ce the tendency tewerd crowdicn migration. Also, 
already mentioned was the fact that incveased crystal size re- 
reduced crowdion migration. Particular attention was paid to the 
proper choice of values for these parameters, in order to cor- 
rectly determine whether or not crowdion migration actually existed. 
This question of crowdion migration was especially critical in this 
simulation, since the binding energy of a particular atom is a di- 
rect function of the nearness of'its neighbors. An atom in a 
split interstitial position feels a more repulsive potential than 
an interstitial that has initiated crowdion migration and 
"stolen' a lattice site and has thus reformed the original perfect 
liAGerce With every atom in a normal lattice site. 

ine ChOolce> or Atom 

Some elements tended to initiate crowdion migration more 


than others. This applied to both choice of lattice atom, and 
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Silotec Of interstitial atom. For instance, crowdion migration 
was much more common in tungsten than in copper. This was due 
to both the size of the tungsten atom and the nature of the crystal. 
Also it was found that different Bienent point defects in the same 
lattice produced varying degrees of crowdion migration. A neon 
interstitial is so small and light that it never initiated crowd- 
10n migration, even when only one layer separated it from the sur- 
face. An argon interstitial initiated crowdion migration at the 
Surface, but not deeper in the lattice. A tungsten interstitial 
always initiated crowdion migration, unless placed at the center 
Grea huge 14 x 14 x 14 enast an Pattee. | this 1S an example of 
increasing crystal size to prevent crowdion migration. In 
general it can be stated: the more massive the interstitial, the 
mOre probable crowdion migration. 
2. Range of Potential 

In surface simulation, the range eerie potential was a 
Gre real factor, Since it determined whether or not an atom could 
"see"! the surface. Because copper has been the standard element 
for lattice simulations, many versions of a copper potential with 
various ranges, have been determined. AS mentioned previously, 
Gw [12] calculated Morse potential for copper that effectively 


had an infinite range (150th nearest neighbor). If this potential 


is truncated at very close ranges, i.e., NNl or NN2, the potential 
is seriously underestimated. This under estimate rapidly dimin- 
ishes as the truncation range increases. Since GW parameters for 


the Morse potential could not be used for NN2 interactions, 
Anderman [14] calculated parameters for a Morse copper potential 


that would approximate GW's results in simulations truncated after 
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NN2. He did this by deepening and broadening the well. Although 
Anderman parameters and GW parameters led to very similar results 
iipam infinite lattice, they led to quite different results in this 
simulation. In general, if the range of a point defect potential 
function overlapped a surface, crowdion migration would take place 
toward that surface, because of an imbalance of forces in the nor- 
mal direction. The effect was a little more complex than this be- 
cause of surface relaxation: if the range of the point defect po- 
tential function overlapped a relaxed surface layer, the slight 
force imbalance would again result in crowdion migration. Accord- 
ing to GGMV, "the machine calculation showed that this atom rapidly 
moved...in a direction determined by minor asymmetries in the 
starting conditions...'". it cis seOpper simulation, an -inter= 
stitial placed in the forth layer with a GW potential range of 
3.1 LU (NN4) caused complete crowdion migration, resulting ina 
copper stub on the surface. An identical run with Anderman para- 
meters for an NN2 potential to a range of 2.4 LU resulted ina 
gO) splat interstitial with minor, damped migration to the sur- 
face. Instead of a stub copper atom as before, four copper atoms 
iamthe surtace layer bulged about .4 LU. 

Another example of a short range potential which demon- 
Strated this lack of ability to initiate crowdion migration was 
the repulsive foreign defect potential, with a range of V3) EU: 
In copper, this potential quickly led to split interstitial 


positions and no crowdion migration for all copper interstitials 


menor Dobe Gcoland, A-N=, Milgram, ., and Vineyard, E.H., 
"Dynamics of Radiation Damage,"The Physical Review, V. 120, No. 4, 
wel23/7, Nov 15, 1960. 


28 





except those placed in the first two layers. In tungsten, even 
this short range potential could not retard crowdion migration in 
mmeetO x LO x 10 lattice. Only in the center of a 14 x 14 x 14 was 
a tungsten split interstitial stable. This stability applied only 
to the short range repulsive potential: a repeat run using the 
standard composite potential with a range of 3.4 LU initiated 
crowdion migration. This increased range enabled the interstitial 
to find minor asymmetries in even a 14 x 14 x 14 lattice. 
3. Energy Damping 

Energy damping was accomplished in this simulation by 
reducing each atom's velocity at the end of each timestep. Two 
methods were eee at first, each velocity component of every 
atom in the crystal was zeroed at the end of each timestep. Later, 
the halving of each velocity component at the end of each time- 
step was employed to save computer time. Ina tungsten lattice 
ee results of both methods were the same; all final positions and 
binding energies were identical. Neither method prevented crowd- 
ion migration. In copper, these two methods led to slightly’ dif- 
ferent results. Although the final position and binding energy 
of an interstitial was AGES identical, and although crowdion 
migration was initiated in both cases (GW's parameters and a 
3.1 LU range were used), the zeroed velocity method had damped the 
migration significantly by the time it reached the surface, where- 
as the halved velocity method caused a complete, undamped mi- 


Gaacion to the surface. 
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een) INGERSTITIAL IMPLANTATION 

All interstitials were placed in the obvious holes ina hard 
sphere, close-packed lattice model. (See Figure 7.) Every 
"hole" in the tungsten lattice had exactly the same geometry; 
1.e., two neighbors 1 LU away; four neighbors 2 LU away, four 
neighbors 3 LU away, etc. The only factors which differentiated 
between these identical holes and thus led to different binding 
energies were: layer number, or lattice depth, and open channel 
direction. An interstitial in the third layer was more tightly 
bound than one in layer two, etc. Also, an atom in a given layer 
could be placed in two ae of holes: one in which the inter- 
stitial was in the BCC (010) open chanrel direction, in which case 
Beieminterstitial could “see” the surface; and one in which the 
interstitial was in the BCC (100) or (001) open channel, in which 
the interstitial could not "see" the surface. Note, however, that 
1f an atom could not "see" the surface, ee there was no difference 
between these two positions, since both have two neighbors 1] LU 
away, four neighbors 2 LU away, etc. Note also that even if these 
two sites are identical and possess exactly the same binding 
energies, the difference might still show up in diffusion pro- 
bababalities: an interstitial in a (010) open channel in the 
second layer must only move two lattice units to escape the cry- 
stal. An interstitial in a (100) or (001) open channel must move 
Pee toecither the < or z direction into an open channel, and then 
2 LU to escape; 1.¢., it must move like a knight in chess. This 


extra step might lead to a different diffusion probability. 


30 





Peete TUNGSTEN LAITICE SELF DEFECT 

The tungsten self interstitials and self replacement defects 
were the chosen standard for this analysis, as explained in the 
in the OBJECTIVE. The tungsten defects could be treated as 
lattice atoms and allowed to interact with all other atoms with 
the composite potential (Method 1); or they could be treated as 
foreign defects and allowed to interact with only a repulsive 
potential (Method 2). A total of three different defect positions 
were simulated: an interstitial in a (010) open channel (int A), 


an interstitial in a (100) or (001) open channel (int B), anda 


replacement atom (rep) in a lattice site. (See Figure 7.) 
1. Interstitials 


AS previously mentioned, all tungsten interstitials initi- 
Pee CLOva:01) Wlcgeation when txcated by method 1. An interstitial 
treated by method 2 also initiated crowdion migration unless buried 
imeme center Of an enlarged 14 x 14 x 14 lattice. Because an 
interstitial that has pushed its neighbors away has a lower po- 
tential than one that has not done so, the numerical values for 
binding energy of tungsten interstitials could not serve as a 
true standard for comparison with W-Ne results. Qualitatively, 
however, much could be learned from the W-W energy levels. First, 
1t was expected that all energy levels of a defect found by 
Method 1 would be negative at equilibrium. Values for the compo- 
Site potential can be either negative or positive, but are posi- 
tive only at very small separations. A negative potential energy 
means the atom is bound in the crystal. As shown in Figure 1, on 
the next page, the first two interstitial levels for the W-W 


reaction, Method 1, were at -3.0 eV and -3.1 eV, corresponding to 


oo 


EV. 


BINDING ENERGY, 
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the two interstitials in the first layer (int A and int B). The 
next two levels were at -5.4 eV and -5.9 eV, corresponding to the 
interstitials in the second layer. Interstitials in the third 
layer and deeper had binding energies ranging from -7.4 eV to 
-8.1 eV. The value of -8.1 eV, labeled "©" was obtained from the 
interstitial placed in the center of the seventh layer of a 
Meer x 14 jattice. Since this value, and all other values for 
Method 1 binding energies were reached after crowdion migration, 
they were all expected to be lower than they would be without 
crowdion migration. The only way the true binding energy of a 
tungsten interstitial could have been found would have been to 
find a tungsten crystal Size large enough to contain the crowdion 
Migration of a tungsten interstitial with the long-range composite 
potential. Compute running time made this impossible. 

Also shown in Figure 1 are the binding energies for the 
Method 2 W-W interstitials. Note, first, Poe they were all 
positive. This was again expected, since the potential equation 
for Method 2 is positive over all space. Note, second, that the 
binding energies for interstitials in the first two layers were 
zero. This was because all purely repulsive atoms in these first 
layers escaped the lattice completely. The other positive levels 
shown, were incomplete, because many interstitial positions did not 
possess stable energy levels. The unstable levels oscillated be- 
cause of significant lattice motion, caused by crowdion migration, 
and measured by a short range potential. The range was so short, 
that significant jumps in binding energy occured when an atom 


moved into, or out of the range of the potential. Evidently, 


32 





small vibrations in atoms with an equilibrium distance of about 

Y3 LU from the interstitial, frequently caused crossings of this 
range limit, adding or subtracting eneray from the binding energy 
each time one of them crossed, thus invalidating many of the inter- 
Stitial binding energies. 

The energy levels shown, however, demonstrated the meaning 
Saeampositive “binding energy". The numbers reflected the "amount 
of repulsion" associated with different positions in the lattice. 
The ordering of the levels, i1.e., higher energies for deeper layers 
was expected. An interstitial deeper in the lattice felt more 
repulsion, because it was surrounded by a greater number of re- 
pulsSive neighbors. Again, the level labeled "©" represented an 
Pmerstitial in a 14 x 14 x 14 lattice: but in this case of a 
Method 2 interstitial, crowdion migration did not occur. 

The concept of a positive binding energy may or may not 
be the actual physical situation, but it is still academically 
valuable. Instead of an atom resting near the bottom of a po- 
tential well, as in Method 1, an atom can be "wedged" between the 
repulsive walls of its neighbors. In both cases, the atom is 
taound’ in the lattice. The ordering, spacing, and other cor- 
respondences between the positive and negative levels validate the 
qualitative use of the positive levels. 

2. Replacement Impurities 

A tungsten lattice with a replacement impurity is a per- 
fect tungsten lattice, but Method 1 or Method 2 could be used on 
the atom.in question. For Method 1, a perfect crystal wasallowed 


fomcelax With time, yielding only negligible motion; and the 
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binding energies of the center atom in each layer was recorded 
(see Figure 1). The binding energies of atoms in the first, 
second, and third layers were-5.3 eV,-6.8 eV, and-7.0 eV respec- 
tively. The subsequent reversal of order for levels corresponding 
to deeper layers is a program anomally, caused by the use of a 
fiommbe depth crystal. Runs on larger crystals indicated that the 
Gug@ers would not reverse in an infinite lattice. The "®" level 
at-8.8 ev is the experimentally determined heat of sublimation [26] 
of tungsten. These numerical values for the binding energies 
were valid as standards of comparison for the Ne-W data, since the 
motion and equilibrium positions of the replacement atoms and their 
neighbors Rane noanly identical, and usvally less ener 23 LU in 
both cases. Note that the binding energies of the Method 1 re- 
placement atoms were lower than the interstitial atom levels. As 
previously stated, this was to be expected, since a replacement 
atom rests in the bottom of a periodic potential well in the 
tierce, whereas an interstitial rests in a higher well, because 
it is nearer to its neighbors than the normal equilibrium sepa- 
ration. Also note that if the entire Method 1 spectrum of binding 
energy levelS were used as s standard of comparison for Ne-W levels, 
then the W-W interstitial levels -should be higher with respect to 
the replacement levels, than shown in Figure 1, because of crowdion 
Magna t1On, 

The Method 2 replacement level labeled "©" corresponded 
to a replacement atom in the center of the fifth layer ina 10 x 
10 x 10 lattice. Note that it was not above the interstitial 


levels, as it should have been by comparison with Method 1]. This 
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was the major weakness of Method 2: because it was a measure of 
repulsion only, and neglected the potential well, it under esti- 
mated the values of the binding energies of atoms whose normal 


position was in that well; i.e., replacement defects. 


E. THE NEON DEFECT IN TUNGSTEN 

The neon atom defect was again placed in any one of the three 
[werrece Positions, labeled "int A", “int B", and “rep'', The neon 
defect could be treated by Method 2 only. The neon energy levels 
are shown in Figure 2, on the next page. Note that the energies 
have been multiplied by a mass correction factor. (See Appendix D.) 

1. Interstitials 

The neon interstitial never initiated crowdion migration, 

but as in the W-W case, atoms placed in the first two layers es- 
caped the crystal, and therefore had zero binding energy. Again, 
the ordering of the levels was a measure of the replusion on the 
interstitial, which increased as the interstitial was placed 
deeper in the lattice. Again the "@©' Jevel was the result of 
Placing a neon interstitial in the center of a 14 x 14 x 14 tung- 
eren lattice. 

2. Replacement Impurities 

Again note that the replacement levels are lower than 

would be expected by comparison to the W-W Method 1 standard. It 
is hypothisized that these replacement levels should be higher 
than the interstitial levels, by comparison to the standard. This 
assumption is valid if a Ne-W potential well exists. It is fea- 
Sible that since neon is almost incapable of binding, that no Ne-W 


well exists. On the other hand, the ionization state of both neon 
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and tungsten in the lattice is unknown, so the existence of a 


shallow potential well is quite possible. 


F. CORRELATION WITH EXPERIMENT 
1. Scaling the Levels and Peaks 

Since no direct way of transforming the Ne-W energy levels 
into correctly scaled negative values exists, an arbitrary linear 
scaling factor between KS's data and the Ne-W levels has been 
used. Note that every level or group of levels corresponded to 
an experimental peak in Figure 2, except: in two places: first, 
the broad peak at about 2000° K had no energy level counterpart, 
but could be assumed to correspond to the closely ordered replace- 
ment levels, shifted above the interstitial levels by the above 
hypothesis. Second, the narrow peak at 4 50°K was without an 
energy level counterpart. Note that three interstitiai ieveis 
had zero simulated binding energy because they had escaped the 
crystal. If, however, the assumption that a Ne-W well exists was 
Maer then Some or alll of these three interstitial locations, 
(i.e., two surface layer positions and the open channel position 
in the second layer) would be stable, bound positions, and would 
be expected to generate an energy level in the vacinity of the 
450° kK peak. If the arbitrary eine of the peaks and levels was 
correct, then the peak at 450° K is proof that a Ne-W well exists, 
since a purely repulsive potential would not allow a first layer 
or second layer open channel interstitial energy level to exist. 
The existence of this well, then, would in turn substantiate the 
shift of the replacement levels above the interstitial levels, to 


correspond to the 2000° K peak. 
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2. Probes of Potential Wells 

Various attempts to substantiate the scaling between the 
levels and peaks were made. No approximation to a Ne-W potential 
well could be justified, and the correspondence between W-W and 
Ne-W results was not complete enough to invert and scale the po- 
tential energy levels to realistic, negative, binding energies. 

Attempts were also made to match both Method 1 and Method 2 
W-W energy levels to KS's data. The Method 2 levels were too 
incomplete, and the Method 1 levels required an unknown arbitrary 
reduction of interstitial energy levels to compensate for crowdion 
migration. Too many alternate reductions were possible to choose 
omenas the correct rescaling. 

Another approach that yielded little information was a 
plot of the differences betWeen interstitial and replacement 
energy levels for each layer in the lattice. 

One valuable method of investigating the nature of possible 
Ne-W negative binding energies was to probe the perfect lattice 
Tieoean interstitial at various initial positions and plot the re- 
fwitanmt potential energy of an interstitial vs. position. Since 
the potential was always positive, the results of this investigation 
were potential wells above the x-axis. Although the positive lo- 
cation of these wells was not realistic, the relative depth of the 
wells was significant. The average depth of a neon interstitial 
well, deep in the lattice, was about 4.2 eV (see Figure 8). The 
graph was made by placing interstitials in positions in the (010) 
open channel, and thus represents the barriers that the interstitial 
must penetrate as it escapes the crystal. Note that the wells wece 


not at the obvious holes in the BCC lattice, but were between layers. 
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In the actual simulation, the interstitial rarely fell into this 
well, but instead pushed its two nearest neighbors away and made 
the initial position the low potential position. Here again, sur- 
face effects and relaxation reduced the tendency for interstitials 
to relax into expected infinite lattice equilibrium positions. 
Slight differences in the final equilibrium positions were not of 
Significant importance to binding energies. The actual numerical 
values for the depths of these positive wells were not necessarily 
scaled properly since they ignored the actual Ne-W potential well, 
and were found from perfect lattice probes at time zero, before 
relaxation. Nevertheless, the depth of a well deep in the lattice 
of 4.2 eV agreed well with KS's prediction of 4.5 ev [8] for the 
desorbtion energy corresponding to ee Oe at the front edge of 
the highest peak. 20 K closely approximated the temperature 
that the arbitrary scaling had assigned to a deep interstitial. 
Note that the initial position of a replacement atom was 

3 Lu from its neighbors, and thus because of potential erosion, 
Memiiacdeal potential of zero eV initially. To climb out of this 
well, about 17.5 eV must be supplied. Although this number was 
inaccurate for the same reasons listed above, it did demonstrate 
that even a purely repulsive potential can predict a greater 


binding energy for replacement atoms than for interstitial atoms. 
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VeecoNCLUSIONS 


An arbitrary scaling has been used to correlate the simulation 
results with experimental data. Although the method was not ana- 
lytically sound, no other avenues of approach to the problem could 
be found that could further justify our hypothesis. 

Satisfaction can be gained, however, from the fact that these 
results compare favorably with known data at many interfaces. Our 
model was a tried and proven one, with many successful sputtering, 
channelling, and similar simulations to its credit. This present 
model invariably behaved ina physically valid manner or a manner 
which could be made physically acceptable by varying the con- 
trolling parameters in the program. Specifically, many previous 
experimental and simulated results for infinite crystals were 
reproduced when simulation took place deep in a large lattice, 
ememeas the (110) split interstitial position for BCC structures. 
The Method 1 replacement levels, if found for a much deeper cry- 
stal, would have asymptotically approached very close to the 
8,8 eV heat of sublimation. The simulated depth of the positive 
potential well for interstitials of about 4.2 ev pieces approx1- 
mated KS's prediction of 4.5 eV. 

All avenues in additional computer simulation have not been 
exhausted. Future simulation of Argon, Krypton, and Xenon defects 
in tungsten should be fruitful. Comparisons between the relative 
locations of these new energy levels might further substantiate 
this Se In particular, if simulation can explain why KS's 
neon data contains five desorbtion peaks, while their Argon, 


Krypton, and Xenon data contain four peaks, it will be a major 
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success. KS also gathered data for different crystal surfaces 
and different angles of incidence, which might be investigated 
by computer simulation. 

This simulation was also important in that it investigated 
the lattice surface; a topic which has not received as much 
attention as infinite crystal dynamics. Since radiation damage 
theory and modern transistor theory 1S very much concerned with 
the crystal surface, the computer simulation field will undoubtably 
increase their emphasis on surface effects, with considerable 
attention toward better ways of treating a foreign interstitial. 

A new exhaustive book which reports on the present state of 
knowledge in all these areas, with emphasis on rere menial re- 
sults has just been published. It iS a report on the proceedings 
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APPENDIASA: CRYSTAL GEOMETRY 


The computer program can call any one of nine lattice generator 
subroutines: three face-centered cubic subroutines (LL100], eran 
merdi)), three body-centered cubic subroutines (C B100], [B110], 
[B111]), and three diamond subroutines ({p100], {p110], (p111]). 
The diamond subroutines are never used and therefore not compiled; 
but provision has been made for their future inclusion in the pro- 
gram. The dimensions of the lattice chosen were controlled by the 
input data variables ae Cry], and (1z]. Each atom in the cry- 
stal was numbered, in the order x followed by z, followed by y. 
Pemeeene surtace layer (Y = 0) of the tungsten 10 x 10 x 10 lattice, 
atoms were numbered from 2-26; for the first layer below the sur- 
face, atoms were numbered 27-51, etc. Atom number 1 was the pri- 
Mary, O1 point defect atom. ierlewesrtnermunver of the Yact mosite 
atom, or the last atom in the eighth layer, number 200; and Paes 
was the number of the last atom in the crystal, number 250. 

The placement of point defects was accomplished as follows: 
after the desired perfect lattice was built to the desired size, 
subroutine PLACE was called. Three types of defects were allowed: 
vacancies, interstitials, and replacement impurities. The type 
and location of the defect were controlled by input data vari- 
ables: the type of defect CiTyPe), an atom number [NVAC] and a 
displacement vector [D1X, Dly, Dlz] in Lu. If [ITYPE] = 1, a 
vacancy was created in site number Cnvac]. This "removal" was 
accomplished by setting (LCuT (NVAC) J = 1 which “turned off" the 
atom, removing it from all calculations. If LITYPE] = 2, an 


. e = ee viep t T rr 
interstitial was created in a position L-DIX, -DIY, + DIZ} LU frem 
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site number [NvAC]. This interstitial was always atom number 1. 
Since atom number 2 was always at the origin, number 1] could be 
placed using a displacement vector from the origin (from [nvac] = 2) 
or uSing a displacement vector from a site next to the interstitial. 
If [ ITYPE] = 3, aA vacancy was created in site number [ nvac] and 

a replacement impurity, put in its place. Note that for both 
fener | = 2 and 3, elther a foreign or self defect could be 

Placed. For the case of either the self-interstitial or the self- 
replacement atom (giving us back the perfect crystal), either 
method 1] or method 2 of calculating the potential could be used. 
The choice of methods was also an input parameter: for method 1, 


C1Q] =i and for method 2, Cra] = 2. 
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Pore Nolxs. POTENTIALS 


(In this appendix and all subsequent appendices, the brackets 
denoting program language are dropped. Program language is still 
written in all capital letters.) 

A. BORN-MAYER REPULSIVE POTENTIAL 
1. Potential Energy: For the lattice atom interactions, the 


Born Mayer potential equation is: 


Vag meno eas Se (1) 


or 


POT = EXP(EXA + EXB*DIST) (1A) 


For bullet-lattice atom interactions, 


vee = exp(A! + nae - ce) (3) 


where aa is subtracted to retard the potential so that it 
goes to zero at the nearest neighbor distance. itn the progiam, 
the V.. equation is: 
oe) 
POT = PXP(PEXA + PEXB*=DIST) - PPTC (3A) 


peemiborce: For the lattice atom interaction, 
~OV.. 
Force = ae = -B exp(A + Br...) 
ae 1j 


se [ (én -B + A) + Br, J (4) 


in the program, &#(-B + A) = (ALOGL-EXB*CVED]+ EXA) = FXA, 


where CVED is a conversion factor for units, and 
FORCE = LFXA + EXB*DIST] (4A) 


For bullet-lattice atom forces, 


-OV.. 
Force = ae = - Biexp(A'+B'r,.) = exp[in-B'+A')+B'r 5 5] (5) 
ij J 
where (In-B'+A')=(ALOGL -PEXB*CVED]+PEXA) = PFXA, and 


FORCE = EXP({PFXA + PEXB*DIST] (5A) 
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For bullet-lattice atom forces for which lattice atoms enter the 
Tange of the Bullet's force during the timestep, a first appr oxj- 


mation to the force is Qliven by: 


Vi 57V, .(ROB) 
ROree = Se (6) 


ree 
15 
(See Figure 9 ). Vig7 V4 jROE) 1s the retarded Potential given 


above. a Poewersian) factor 1S needed to Preserve the Proper units 


meet On (3A)- (PEXA) becomes (2n CVED+PEXA) = PAC 
aa) EXP[PAC+PEXB*pIsr] and V;j(ROE) = perc becomes PFPTC, which 
1s V5 j(ROE) calculated after the Pac Substitution. Finally, 


Tiq7ROE = DIStT-Rop = DFF, and 


1) 
( +PEXRB* = 
FORCE = EXP[ pac oe DEC PFPTC (6A) 


potem( A) is useg PereO= DIST Ss pon _ DTI = ROEM, and 


fo) as used for ROEM < prsqt < ROE, 


* 
v 
o.° Pe 
s 


B. MORSE POTENTIAL 


1 Potential Energy: LOD Wateics atom interactions, the Morse 
. = , SP es 


POtential equations is. 


oe OE St Boel ale oe (2) 


tl 


expl (2np + OO 9) (20)r 5 J~expl (on 20) 4 Wo) (A)es 


= Expl ( ALOG( Deon) + 2.*ALPHA*RB) AS *ALEHA*CVR) *DIST) “ee 
- EXPL (ALOG(2.*Dcon) + ALPHA*RE) ~(ALPHA*CvR) *prgq] - a 


ExP( OGD1-cGB1*prs7]- EXPL CGp2-ceB2*prs7} | (2K) 


2. Force. POR tat tice atom interactions, : 


-09.. | 
Hor Gay = aes = pL 2a empl -20(r 5 ry} 29 expt -O(r 5 sry” 
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a) 





= expL0n( 20) + (@nD+20x -(20)x 5 J-explona +(&n( 2D) +r ) -(Q)r i od 


= ExP({ ALOG( 2*ALPHA*CVR*CVED) +ALOG( DCON) +2. *ALPHA*RE 
-(2.*ALPHA*CVR) *DIST] 

- BEXP[ ALOG( ALPHA*CVR*CVED) +ALOG(2.*BCON) +ALPHA*RE 
-( ALPHA*CVR) *DIST] 


EeEXPLALOG( -CGB1*CVED) +CGD1 ) -CGB1*DIST] 


J 


- EXPL (ALOG( -CGB2*CVED) +CGD2 ) -CGB2* DIST] 


=m CGhil-CGRL*DIST |= EXPl CGF2-CGB2*DIST | 


wae CUBIC FIT 


(7A) 


1. Potential Energy: The best cubic fit between the BM and 


Morse potentials is calculated in Subroutine CROSYM. 
tential equation, defined between ROEA and ROEB, is 


3 Z E es 
POme- “Coote, FteCPZtee. + CPbr. . +°CPy 
1) 1) 1 J 


or POT + DIST*(DIST*( DIST*CP 3+CP2)+CP1)+CP@ 
-OPOT 


Ones 
AJ 


2 
2. Force: Force = = =3CP3r -. =2CP2r ...<=-CP1 
a) 1) 


i 


(-3.*CP3*CVED) x, +(-2.*CP2*CVED)x ; 5+ (-CP1*CVED) 
2 
Cite eh lites. ONY, OF 
1) 1) 


FORCE = DIST*(DIST*CF2+CF1)+CFG. 
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(9) 
(8A) 
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PPP eNDIX C: AVERAGE FORCE METHOD AND TIME DURATION THEORY 


A. AVERAGE FORCE METHOD: The average force technique has been 
explained in great detail in Ref. 15. 1t was summarized in 
Chapter 3, and therefore discussion here is limited to the average 
force method in the program language. 

Paemeene desired lattice 1s built, the position of the ith atom is 
eeeued Simultaneously in RX(1), RY(1I), RZ(1); RXK(1), RYK(I), 
RZK(I); and RXI(I), RYI(I), RZI(1I). The latter set of coordinates 
never change and are used for comparing new positions to original 
Peattions, and in calculating DX(I), DY({I) and DZ(1I) for output. 
The middle set: of coordinates containing the letter K are for 
storing the initial positions at the beginning of each timestep. 

A step by step summary of the average force method, showing X 
coordinate calculations only, roilows: 

imgeebased On the position, RX(1I), of the ith particle at the be- 
Oinning of the timestep, the force FX({I) is calculated in STEP. 
Eee ( E) 1s Stored in RXK(I). 

3. The new, temporary position, RX(I) is calculated, based on 
mhemrorce at RXK(1): 


yt 


t 


2 , 
Were vAtet (At) 72em (10) 


Or 


RX(I) = RX(1)+DTOD*( HDTOM*FX(1I)+VX(I). (10A) 

4. A new force is calculated, based on this new position RX(I). 
5. VSS stores VX(I), the original velocity of the ith atom. This 
velocity is half the velocity of the ith atom in the previous 


1, factor being an arbitrary damping multiplier. A 


timestep: the 
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| 


new velocity, based on the new force, is found: 


y 


V + FAt/2m (ala 


OT 


VX(I) = VSS+HDIOM*FX(I). (11A) 
6. The final position is calculated, based on the average of 
these two velocities 


xt 


xt bat(V+V") (12) 


or 


N 


O10) RXK (I) +(VX(1I)+VSS ) *HDTOD (12A) 
The resultant velocities are halved, a new timestep duration is 


calculated, and the process repeated. 


9 


Bae liMESTEP DURATION THEORY 
This simulation uses the best possible estimate of a timestep 


fumeition, PT, as calculated from the PLCcene sta te OmminesLOLces 


) 


and energies in the lattice for use in the next timestep. To 





limit motion to an increment small enough to preserve the accuracy 
of the average force approximation, we define DII as the maximum 
distance any atom is allowed to move in one timestep. 


From (10), we find 


AX. = (V.+ F.At/2m)At. 
ot ab - 


Therefore, 


At AX. /(V,+F,At/2m) . (13) 


For 


V, 7? F.dt/2m, At = AX,/V,. (14) 


If we find the fastest moving atom and assure that it does not 
move more than DI, we have limited the motion of all other atoms 
to less than DTI. 
Thus, 


DT = (DTI*CVD)/EMAX = FDIT/EMAX (144) 
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where 
EMAX = SQRT(VX(1I)*VS(1I)+VY(I) *VY(I)+Vz(1I)*Vz(I). 
For V. 27> F.At/2n, 
1 i 
emAXx. 
At = AX./F ,At/2m ae 5 (15) 
1 
Anatagous to above, we find the most stressed atom and assure 


that it does not move more than DII. Thus, 


DT = sSQRT( (2.*PTMAS*DTI*CVD) /FMAX] 


SORTL TFAC/FMAX] (15A) 


I 


where FMAX= SQRTLFX(I)*FX(I)+FY(I)*FY(1)+FZ(1I)*FZ(I)]. 

Since rigorously we cannot make either or these limiting assumpt- 
Benemewe must GO back to Our Original equation for DI, equation (13). 
Since this equation involves DT, we proceed as follows: 

1. Assume Vi << F,At/2m SiidmeartenlLeate At irom (415A). 

2. Insert this preliminary value for DT in (13) and compare 
Vv, to F.At/2m. If V, is larger, calculate DT from (14A). If 
F,At/2m is larger calculate DT from (15A). 

A complication arises when a foreign impurity is in the lattice, 
Beeause Of a4 Variation in the value for m in (15). This is 
especially acute when the differences in masses are great. The 
method used to solve this problem is as follows: 

1. If either FMAx = F 


OEP EM = Vv the entire proceedure, above, 


1 1’ 


is followed using the mass of the bullet for m. 

Ze it) beth FMAX 7 F, and EMAX 7 V,, the entire proceedure is fol- 
lowed uSing the mass of a lattice atom. 

The ae aaenene that the bullet mass be used if either EMAX or 


FMAX describe the bullet circumvents the problem of having the 
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bullet the fastest moving atom, but not the most stressed aton, 
or visSa versa. 
., . -14 , 
To begin the problem, an arbitrary value of 10 seconds is 
assigned to DI. If at any time in the program EMAX = FMAX = zero, 


10714 1S again assigned to DI to prevent division by zero. 
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PrueNDix DD: SUBR@UIINES STEP, ENERGY, AND LOCAL 


Pee DLSLANCE CALCULATIONS 

In all three Simecoeinee Sle wee Ne GY wand LOCAL, a method of 
finding all atoms within a given radius of another atom was needed. 
For lattice atom interactions, atoms inside ROEA, ROEB, and ROEC 
were found; for the foreign interstitial interactions, atoms inside 
ROE were found; and for LOCAL, atoms inside ROEL of a point defect 
were found. The time Saving technique used to do this was to 
successively eliminate all atoms with an x component difference 
greater than the given radius, then similarly for y components, 
then for z components. The resulting volume not eliminated is a 
cube circumscribing the desired sphere. Finally the time-con- 
suming test of eliminating all atoms for which the desired radius 


-_ RP. 


is less than SQRT (DRX*DRX + DRY*DRY + DRZ“DRZ) is applie 


(2. 


to 


only atoms inside the cube. 


B. SUMMATION INDICES IP AND IQ 
Interactions We are and cer are found by evaluating all values 


in the half matrix. For example F are found, then 


Meats a Te 
Fi 3? Pay? ne tC. ie variable IP controls the starting point 
for the j summation. IP is always set to I + 1 to avoid the re- 
petition of finding ae and ce T@ controls the starting pont 
for the i summation. If the primary is to be treated as a lattice 


atom. ©FO = 1. If it is to be treated as a foreign particle, 


1O>—=2,,and all Pa are found separately. 


Cae DOU STRIBULION OF FORCES AND POTENTIAL ENERGIES 


The forces Pa are equal and opposite on i and j: 1.e., F.=~F 
a 


The potential energies are split in proportion to the reduced 


ae 





mass of the interacting particles. This 1s easily understood by 
Observing the Kinetic energy distribution of an elastic collision 


of mand M where M7? m. We find that m carries away almost all 


{M 


() ETor Hror’ 





the kinetic energy: specifically it carries away 
If a pair of atoms are to behave elastically, the potential 
energies which are transformed into kinetic energies of motion 


must be split in the same manner. For this reason, 


TMAS 
=a ee 3 am ¥ e 
PPE(1) (5 oe POT = BSAVE*POT; 


BMAS 
EEE) ae POT = TSAVE*POT. 


note, for BMAS = TMAS, BSAVE = TSAVE = 4, and the energies are 


Split equally. 


Dee oU OK OULINE LOCAL 

LOCAL measures the change in potential energy associated with a 
Sphere of radius ROEL surrounding a point defect. It sums up the 
potential energies of each atom found inside thissphere at time 
zero. It remembers these atoms, and for each timestep re-sums 
the potential energies of these same atoms. The sum, total local 
potential energy TLPE, is subtracted from TLPE at time zcro 
(TLPEZ), to give a measure of the change in potential energy 


(DLPE) inside ROEL. 
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APPENDIX E: COMPUTER PROGRAM GLOSSARY 


fee. in this glossary, the terms “point defect atom", "bullet", 

and "primary" are synonymous; and the terms "lattice atom" and 

"target" are Synonymous. 

ALPHA: Input Morse potential parameter 

BSAVE: Target mass/(target mass + bullet mass); distributes 
potential energy between target and bullet 

BIND: Negative of the total potential energy (TPOT) at time zero 

BMAS: Mass of bullet in amu 

Bumeetl: Alpha=numeric array for point defect material 

CFO, CF1, CF2:- Force parameters of cubic fit between Morse and 
Born-Mayer functions 

CGBl1, CGB2: Morse potential parameters 

CGD1, CGb2: Morse potentiai parameters 

CGFl1, CGF2: Morse force parameters 

ieeerecr i, CP2, CP3: Potential parameters of cubic fit between 
Morse and Born-Mayer functions 

Gyp; CVR x tone converts lattice units to meters 

GVE> 1.6 x tou, converts electron volts to joules 

CVED: CVE/CVD, a ratio used to avoid repeated division 

CVM: 1.672 x Ome converts Ae oiiic mass units to kilograms 


CVR: LU in angstroms; converts lattice units to angstrom units 


Piece DIY, DIZ: Displacement coordinates for location of interstitial 


from reference atom, NVAC 
DCON : Input Morse potential parameter 
DFF: ROE-DIST, the distance closer than ROE that an atom is to 


the primary 
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DIST: 


PEPE: 


DTOD: 


DIOM: 


DIOMB: 


Distance between any two atoms 
TLPE-TLPE%, the change in total local potential energy 


Simce time zero 


Die ORZ: X,y¥,2 Components of DIST 


Length of a timestep in seconds 

Number of lattice units most energetic atom may move in 
one timestep 

DT/CVD--a ratio used to avoid repeated division 
DT/PTMAS--a ratio used to avoid repeated division 


DT/PEMAS~--a ratio used to avoid repeated division 


Bem, PY(1), DZ(1): Change in position of ith atom from initial 


MAX: 


EV: 


EVR: 


position at time zero 
The maximum energy encountered in any cycle 
Ieee meta si c.cC.r on, VOlts 


Primary energy in kilo-electron volts 


EXA, EXB: Input Born-Mayer potential function parameters for the 


F2: 


FA: 


FDTTI: 


FM2: 


FMAX: 


FOD: 


target 
Square of the force on a specific atom 


The component force increment on an atom 


DTI X CVD, a parameter used to determine DI my maximum 


energy method 

A small number used in checking potential energy zero 
point 

FM squared 

Maximum total force on the most stressed atom in the 
crystal 


FORCE/DIST--a ratio used to avoid repeated division 
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FORCE: 


Numerical value of the force function with a variable 


parameter 


feet), FY(I), FZ(1): x,y,2 components of total force on an atom 


rxXA : 


HBMAS : 


HDTOD: 


HDTOM: 


HDTOMB: 


HTMAS: 
ae: 
ss 
IDEEP : 


THI: 


1 
26) 
ho 


ES : 
as: 
LH ; 
WIVAY : 
iN: 


iP . 


10; 


ESHUT: 


ri 


TIT: 


Born-Mayer force function parameter 


*3 BMAS--a ratio used to avoid repeated division 
44 DTOD--a ratio used to avoid repeated division 
*s DTOM--a ratio used to avoid repeated division 
*4 DTOMB--a ratio used to avoid repeated division 


*% TMAS--a ratio used to avoid repeated division 


Variable in cubic fit subroutine 


Number of mobile layers 
Alpha numeric array for program title 

Morée functicom per ameter s 

bullet element 

Pvoc and Oldentati1on of crystal 
target element 

Same as IDEEP 

Odd-even integer used to determine atom site establishment 
Subscript value of atom. Used in subroutines STEP and 
ENERGY 

or not a self defect is 


Parameter that determines whether 


to be given a repulsive potential or a composite attractive- 
repulsive potential 

A parameter used to shut down the program 
Unscaled fixed point x coordinate used in lattice generation 


Odd-even integer used to determine atom site estab.iishr.ent 
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TYPE: 


a, IY, 
wie: 
oe): 
ot: 
oS; 
oar: 
KF; 


KT: 


Pevr( 1): 


LD: 


| 


Parameter used to determine the type of point defect: 
vacancy, interstitial, or replacement 

PZ euUnper Orc yz planes of crystal 

Variable in the cubic fit subroutine 

Parameter in the BCC(111) lattice generation subroutine 
Unscaled y coordinate used in crystal generation 


Variable suSsed to eStablish atom Sites 


Final K in LOCAT (K) assigned to an atom 
Unscaled z coordinate used to establish atom site 
Used to identify an ith atom which is not included in 
calculations | 
The highest numbered atom in the mobile layers 


PiewiigieSe MuUMDCI ea ATOM 1n the entare crystal 


LOCAT(K): Dimensioned variable that remembers the numbers of the 


LS: 


MCRO: 


ND: 


NEW: 


PEeaGLe - 


NRUN: 


atoms within a radius ROEL of the primary at time zero 
Variable associated with each of the nine lattice 
generator subroutines 

One number higher than the order of the fit between the 
Born-Mayer and Morse potentials, always 4 in this simu- 
lation 

Data output increment, in numbers of timesteps 
Parameter used to determine whether or not atom numbers 
have been stored in LOCAT(K) 

Page numbering variable 

Parameter used to determine whether or not to read 


additional data cards 
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NS: 


ING: 


Nil: 


NVAC: 


EAC: 


PBMAS: 


PEXA, PEXB: 


PPeic: 


PFXA: 


PKE(I): 


eeANE : 


POT: 


PPE(I): 


PriIc: 


PTE(I): 


EMA: 


RE: 


RO: 


ROE: 


ROE2: 


ROEA : 


ROEB: 


ROE GC: 


ROEC.- 


ROEL : 


Initial print statement timestep number 

Timestep number 

Timestep number limit before shutdown 

An atom number used to establish point defects or used as 
a reference point for interstitial placement 

Bomameter for bullet) force function correction 

Primary mass in kilograms 

Input Born-Mayer potential function parameters for 
the bullet-target interaction 

Primary force function evaluated at ROE 

Primary force function parameter 

Kinetic energy of the ith atom 

Alpha-numeric array for lattice orientation 

Sect veoeme. 2 tons 
Potential energy of the ith atom 

Primary potential function evaluated at ROE 

Total energy of the ith atom (potential + kinetic) 

Target mass in kilograms 

Input Morse potential parameter 

Spacing constant in FCC(110) lattice generation subroutine 
Nearest neighbor distance 

ROE squared 

Maximum cut off for Born-Mayer potential 

Minimum cut off for Morse potential 

Maximum cut off for Morse potential 

ROEC squared 


Radius inside of which local potential energy 1s iund 
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ReEte: ROEL squared 


ROEM: ReOb-PET, Legion in which modification of repulsive force 
must be made 

RX(I), RY(I), RZ(1I): x,y,z coordinates of an ith atom at any time 
Pee), RYI(1), RZI(1): x,y,z coordinates of an ith atom's initial 


position 


Peet), RYK(1), RZK(1): x,y,z coordinates of temporary position of 


SAVE: 


ee. SCY, SCZ: 


an ith atom during force cycle 


= POL 


Mv COOrGdinatesscale factors 


SSeZ: nmemocale faclLor used for the FEC(111) lattice generator 
subroutine 

START: An optional timing variable, not used in this simulation 

SUM Waraablo wm cubic £11 subsmenmdane 

TARGET: Alpha-numeric array for target material 

TSAVE: Bullet mass/(target mass + bullet mass); distributes 
potential energy between target and bullet 

ve: Total energy of all crystal atoms (kinetic + potential) 

TEMP : Temperature of lattice in degrees Kelvin. Not used in 
this simulation 

TFAC: A time factor ratio used to determine DI by maximum force 
Bead 

TFACB: TFAC for the bullet 

THERM: Thermal energy of atom. Not used in this simulation 

TIME: Elapsed problem time in seconds 

DEP ES: Total local potential energy of atoms within a radius ROEL 

TLPE@: TLPE at time zero 
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TMAS: 


TPKE: 


TPOT: 


VSS: 


Target atom mass in amu 


Total kinetic energy of all crystal atoms 


Total potential energy of all 


Storage variable for velocity 


Pee ye vY(l), VZ(L): x,y,2 components 


X, Y, Z: Unscaled coordinates used in 


crystal atoms 
components 
of ith atoms velocity 


crystal generation 


faerrcet jy: Relaxation in -y direction of ith layer in L.U. 


Ze: 


Hioatang point Lorm of JIT 
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FIGUPE 7. 


BCC (100) ORIENTATION 
HARD SPHERE MODEL 


SOLID LINE: ‘°“Y° PLANE 
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FIGURE 8. 
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FIGURE 9. 
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